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Abstract 

The energy loss of charged particles in silicon can be approximated by a simple analytical model. With help of 
measured charge deposits in individual channels of hit clusters their position and energy can be estimated. Deposits 
below threshold and saturated values are treated properly, resulting in a wider dynamic range. The proposed method 
gives improvements on both hit position and energy residuals. The model is successfully applied to track differential 
energy loss estimation and to detector gain calibration tasks. 
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1. Introduction 

The identification of charged particles is essential in several fields of particle and nuclear physics: particle spectra, 
correlations, selection of daughters of resonance decays and for reducing the background of rare physics processes 
|[l[|2l. Silicon detectors can be employed for identification by proper use of energy deposit measurements along the 
trajectory of the particle. The aim of this study is to provide a simple method to evaluate energy deposits and other 
derived quantities for a wide range of particle momenta (Jiy - p/m - 0.56 - 10.0) and detector thickness, based on 
the precise knowledge of the known underlying physics processes. 

This article is organized as follows: the model is motivated and discussed in detail in Sec.|2] Among the various 
applications we will deal with the estimation of hit position (Sec. O, hit energy deposit (Sec. |4]i, differential energy 
loss for tracks (Sec.|5]l, as well as detector gain calibration for tracks (Sec.|6]l. This work ends with conclusions. 

2. Simple energy loss model 

When a charged particle traverses material it loses energy in several discrete steps, dominantly by resonance 
excitations ((J-function) and Coulomb excitations (truncated term). This latter is the reason for the long tail 
observed in energy deposit distributions. It is generally difficult to describe the process with a simple function. 
The most probable energy loss A, and the full width of the energy loss distribution at half maximum Fa [3J can be 
approximated by 



A = ^ 
Fa =4.018^ 

where 



log 0.2000 -^^-5 



(1) 
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Figure 1 : Probability density functions for particles with j3y = 1 .00, at given path lengths of 20, 50, 100 and 200 jum. 



is the Landau parameter; K - AnNAr^nieC^ - 0.307 075 MeV cm^/mol; Z, A, p and / are the mass number, atomic 
number, density and thickness of the material, respectively; and 6 is the density correction 

The probability of an excitation, energy deposit, along the path of the incoming particle is a function of /3y = p/m 
of the particle and depends on properties of the traversed material. The conditional probability density p(y\l), deposit 
y along a given path length /, can be built using the above mentioned elementary excitations combined with an 
exponential occurrence model. The details of the microscopical simulation for silicon can be found in Refs. ||3|], ^ 
and (@]. The result of these recursive convolutions is a smooth function (Fig.[T]i. 

In the experiment the deposit values y are measured. For a fixed y the conditional probability can be regarded 
as function of /, thus p{y\l) = fy{l). Values of /,,(/) for particles with fSy - 3.16 are plotted in Fig. |2] where the 
values corresponding to y = 50, 100, 200 and 300 keV are shown. Other ySy settings show very similar behavior (J3j - 
0.56, 1 .00 and 10.0, not plotted). Apart from a multipUcative factor, they are well approximated by a combination of 
simple functions, exponential and Gaussian: 

if / < I* 

(2) 




if z > r 



where the limit is 



r - ip - vcTi. 



The peak position Ip and standard deviation cr; are both functions of y and Py, v » 0.65 is constant. Note that fy 
is constructed such that the value and the derivative are continuous at the limit /*. The peak position Ip is roughly 
proportional to the deposit y, their relationship can be approximated as 

y^slp[l+Omiog{lp/lo)] (3) 

where /q is some reference path length (Fig. O. In this work Iq - 300 ;um was chosen. At the same time cr; has an 
approximate first order polynomial dependence on y (Fig.|4]i. 

As it was shown above, some interesting connections between the deposit and path length can be observed, for 
wide range of fiy values. Unfortunately the dependence on ySy, and /, is still present. In the following (Sec. 12.1b we 
will show how that can be suppressed or even eliminated. 
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Figure 2: Fits of the conditional probability density /, (/) as function of the path length / for particles with j}y = 3.16 (see Eq.^2)). The simulated 
values (points) and the fitted curves (lines) con'esponding to y = 50, 100, 200 and 300 keV are shown, from left to right. 
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Figure 3: Relationship between the measured deposited energy y and the position of the path length peak /,,, shown together with fits using Eq.[3] 



2.1. Most probable differential energy loss 

The path length dependence of the most probable energy loss A (Fig.|5]l has a form that is very similar to Eq. (O. 
The dependence can be approximated by 

A(0*eZ[l+fllog(///o)] (4) 

where s is the most probable differential energy loss along a reference length Iq, hence s = A(/o)/Zo- For the relevant 
/3y region studied in this paper a ^ 0.07, to a good approximation^ For some applications, such as hit position 
estimation, the assumption of linearity may be enough (see Eq. (fTTl i later). 

Hence it is worthwhile to check whether the dependence on /3y and Z can be suppressed by taking the most probable 
deposit A as variable of the conditional probability. Values of p{y\A) = /v(A) for particles with /3y - 0.56, 1.00, 3.16 



Based on Eq.^ with some approximations, the expected mag nitude of a would be 1/ logdm^c^ Kz^ ■ Z/A ■ p/I') » 0.064. 
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Figure 4: Relationship between the measured deposited energy y and the standard deviation o"; of / distribution. 



and 10.0 are plotted in Fig.|6] where the values corresponding to y - 50, 100,200, and 300 keV are shown. The 
distributions for a given y value are remarkably similar for all /3y values. 

The conditional probability density for a given y, as a function of A can be described as 

(■exp(^ + ^), ifA<A* 
MA)^C(y)-\ ^tv)n ...... (5) 

where C(y) is a universal function of y. Note that fy is constructed such that the value and the derivative are continuous 
at the Umit A* that is given by 

A*(y)=y-y(TA(y) (6) 

where v ~ 0.65 is constant. Fig. |7]-left shows that cr^ is a first order polynomial of y, practically independent of ySy, in 
the form 

o-A(y) ^o-Q+by (7) 

where ctq ~ 2 keV, b » 0.095. According to Fig. |7]-right the y dependent coefficient of Eq. (|5]l is very well approxi- 
mated as 

C(y)c^l/crdy). (8) 

The parameters of the above detailed energy loss model and their short explanations are listed in Table [T] Note 
that these four numbers are the only parameters of the description. Usually the detector and readout noise can be 
neglected. If this is not the case the term cr\ in Eq. (|5]l and ([8]) should be replaced by cr^ + cr^ where cr„ is the standard 
deviation of the Gaussian noise. 



Table 1 : Parameters of the energy loss model. 
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0.095 


linear coefficient of cr/^(y) 
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Figure 5: Relationship between the most probable energy deposit A and the path length /, shown together with fits using Eq.|4] 



2.2. Left truncation and right censoring 

During readout the deposited energy is converted to measured ADC values through several steps: primary and 
secondary electron-hole pairs, current signals, front-end electronics, digitization |7]. In this study we assume that the 
response of the detector system is linear in the threshold-to-censoring region. Signals below threshold are truncated, 
since they produce no output. Signals above a certain level are censored, hence only the fact that the deposit was 
above that level is known. Note that the censoring level does not necessarily coincides with the saturation level of the 
readout electronics: the former is simply chosen as the limit of linearity. 

If the measured value y is below or above a Umit r, the corresponding values can be calculated by integration 



/v<r(A) 



-L 



fy>,(A) 



-r 



/v(A)y. 



Note that since cr^ is a function of y, the integrals are difficult to perform. The simulated values with several truncation 
thresholds and censoring levels are shown in Fig. [8] They are well approximated by the following functional forms: 



/v<r(A) = 



fy>,(A) 



1, 



if A < r - cta 
if A > f - cta 

if A < f -H (Ta 
if A > r -H cta. 



(9) 

(10) 



2.3. Log-likelihood minimization 

In summary, for a given measured deposit y the probability density will be a function of A only. The measure of 
goodness can be obtained with help of the log-likelihood value as 

For a given channel, with measured {y) or limited deposit values (f), the corresponding;^^^ can be derived from Eqs. (|5]l- 
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Figure 6: Fits of the conditional probability density / (A) as function of the most probable energy loss A for particles with fiy = 0.56, 1.00, 3.16. 
and 10. (see Eq. (5)). The simulated values (points) and the fitted curves (lines) corresponding to y = 50, 100, 200 and 300 keV are shown, from 
left to light. 



(H) and (l9ll-(fT0ll and they are shown below: 



-2v(A-y) 




if A < y - vo-a(3') (11) 
if A > y - vcrA(3;) 



if A < r - crAt) 

^i,(A)=r*^"'- (13) 

■ [0, ifA>f + crA(f) 

where cr/^(y) is given by Eq. The terms on the right side contain only linear or positive definite quadratic functions 
of A. We can estimate e for single hits or for the whole particle trajectory by minimizing the sum of corresponding 
chi-square values. 
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Figure 7: Left: Relationship between the measured deposited energy y and the standard deviation cta of A distribution. Right: The measured value 
of C(y) for particles with /Jy = 0.56, 1.00, 3.16, and 10.0. The points and curves are successively multiplied by factors of 2 for clarity. Note that 
the plotted function I/ctaO') is the same in all cases except for the multiphcative factor mentioned before. 



The proposed model has several applications, such as 

• for a given hit the estimation of position P and energy deposit via optimizing s and the individual path lengths 
li in the various sensitive units, channels, of the detector (Sec.|3]and|4]i. 

• for a given trajectory the estimation of most probable differential energy loss via optimizing s (Sec.|5]l. 

• detector gain (cross-)calibration with tracks via the variation of the chip by chip gain, hence by modifying the y 
values (Sec.|6]i. 

2.4. Estimation of parameter errors 

In case of normal approximations, the inverse of the observed, sample based. Fisher information I can be used to 
estimate the co variance of the fitted parameters V fl: 

1 LiXf 



I 



2 dajdak 



V 



-r 



where YjiX^ denotes the joint chi-square of hit channels or trajectory hits. 

2.5. Generation of random deposits 

If the expected shape of the distribution of the estimator has to be determined using Monte Carlo simulation, the 
generation of deposits should be fast but precise. The energy deposit distribution function 



1 



o-Aiy) 



(A-.v)' 



if y < A + vo-A(y) 
^+^^1, if3'>A + v...(3.) 



does not have a closed antiderivative, hence a simple inverse transform sampling is not possible. Due to its long high 
energy tail the acceptance-rejection method is also highly ineffective. 

Let us suppose that we have generated a random deposit yo at a given Aq. The value yo can be simply transformed 
to get a proper random deposit y for A with a linear transformation as 



y- 



yo ■ 



o-a(A) o-a(Ao) 
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Figure 8: Fits of the conditional probability /,g,(A) for measurements below threshold (left) and above censoring level (right) as function of A. in 
case of several limiting values. Examples thresholds are 5, 9 and 13 keV, while censoring levels are at 80, 135 and 180 keV. The siinulated values 
(points) are shown for particles with jiy = 3.16, together with the fitted curves (lines, from Eqs. (9) and UOt ). 



The distribution ofy can be deduced from P(3'()|Ao) and it can be shown that it equals with P{y\K). The key of the 
proof is the identity 

CrA(Ao) (Ti^iy) _ ^ 
O-a(A) crA(3'o) 



3. Estimation of hit position 

During the passage of a charged particle through sensitive silicon volumes one or more channels (pixels or strips) 
are hit. The charges created in the adjacent channels are recorded (electrons or holes), they provide input for hit cluster 
recognition. Here the task is to estimate the position P of the hit, location of the trajectory in the central plane, using 
the measured deposits. 

Throughout this paper it is required that the local direction of the trajectory, hence the vector of the projected 
passage A (projected onto the surface layer by drift), is known (Fig. |9l). In other words only the re-estimation of 
cluster parameters is attempted. 



3.1. Standard estimation 

The simplest way of hit position estimation (labeled as Weighted) is the weighted mean of channel positions /?;, 
where weights are the corresponding energy deposits y,-: 

■« Weighted — "^^i ■ 

Another widely used standard reconstruction technique (labeled as First-last) deals only with one dimensional 
clusters jgi fioll . This treatment is natural for strip detectors. For pixels the cluster is projected onto both directions {x 
and y) and the two projections are analyzed separately. Only the first and last, projected and summed, channels are 
used for position determination, because that choice reduces the sensitivity to fluctuations in energy deposition. If the 
energy would be lost steadily, the estimated projected hit position would be 

_ Pf + PL yi-yp , 

J^First-last - Z r TT, \^eff 

2 2(yL+yF) 
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where the x and y components of are the sum of the path lengths in the two corresponding edge channels. The 
above formula can be further corrected for the drift direction if the silicon was placed in magnetic field. 

Due to the use of projections, relevant informations corresponding to the shape of hit cluster and its deposits are 
lost. In this work a new method is presented which takes each individual channel of the cluster into account. 

3.2. Estimation using the model 

The passage of a charged particle through the silicon layer is given by the entry and exit points of its trajectory. 
Even in magnetic field or for very slow particles the passage is very well approximated by a straight line, the entry 
and exit points being on different sides of the layer The deposited charge along the passage is projected onto the 
surface layer by drift, as result of electric and magnetic fields. The local direction of the particle at a silicon unit can 
be calculated using the parameters of the trajectory and the vector of the projected passage A is known. 
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Figure 9: Example cluster with measured deposits (rounded numbers in keV) and the the projected passage X (arrow) at some stage of the mini- 
mization. The pixels below threshold are indicated with a lighter fill (below thresh.), while the outlier pixel has a darker fill than the others. 

The following discussion deals with pixel silicon detectors, but it can be easily applied to strips by taking one of 
the pixel dimensions to very large values. During clustering neighboring pixels (those with common edge or vertex) 
are grouped to form a cluster, a reconstructed hit (Fig.|9]l. For the analysis of the corresponding cluster two types of 
pixels are of importance: 

• A touched pixel is touched by the projected passage, regardless of its measured deposit. It can belong to the 
cluster (non-zero deposit) but it may be an empty pixel as well, if the deposited charge was below the threshold. 
The intersections of the projected passage and the edges of a pixel naturally give the projected path length inside 
that pixel. 

• An outlier pixel is not crossed by the projected passage but it is a member of the cluster (non-zero deposit). 
In other words this pixel is an outlier, mostly left by a secondary particle, often 5-electrons, created during 
the passage of the charged particle in the silicon. The projected passage does not have any section inside the 
pixel, but it is possible to construct a measure of the pixel-passage distance, a sort of negative path length (see 
Sec. 13.3b . That notion is important in order to properly take into account the effect of charge diffusioro and to 
be able to move the projected path towards the outlier. This extension also ensures a smooth convergence of the 
minimization. 
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Figure 10: The three possible relative placements of projected paths (line segments) and an outlier pixel (filled rectangle). The distance of the line 
segment and the rectangle is shown by the corresponding arrows (/la,4,c). 




L 



Figure 1 1 : Line segments (dash dotted lines) on the surface of a pixel (filled rectangle). They can be bound by endpoints (open circles) or crossing 
points, a: The length is insensitive to small change of the hit position, b: The length changes when the hit is moved, c: Normal vectors ni, ni and 
the vector of the projected passage X. 

3.3. Projected path length 

The task is to determine the projected path lengths At for each pixel /, for a given hit position P. Using the 
coordinates of the two endpoints, find points where the track crosses the horizontal (along x direction) or vertical 
(along y direction) pixel division lines. Determine the dominant direction of the cluster by choosing longer projection 
of the passage to the axes. Sort the endpoints and crossing points along their coordinate in the dominant direction. 
The lengths of the resulting line segments will give the projected path lengths /I, for touched pixels. In case of outUer 
pixels determine the distance of the projected Une segment to the pixel (Fig.fTOli by evaluating the combinations 

• two endpoints of the projected path line segment with respect to the four vertices and four edges of the pixel 
(cases a and b), 

• projected path line segment with respect to the four vertices of the pixel (case c), 

and choose the negative of the smallest distance. The outlined procedure ensures that by varying the hit position P, 
the projected path lengths of touched and outlier pixels will change continuously, without jumps. 



-In case of 300 fim thickness its contribution is about 5 fim. 
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3.4. Derivatives 



How does the projected path length A change if the hit position P is varied? The corresponding line segment can 
be bound by endpoints or crossing points. The length is insensitive to small changes of the hit position if both points 
are endpoints or if both points are crossing points on opposite sides of the pixel (Fig. [TT]a). The length will change 
only if one of the points is an endpoint or if both points are crossing points on neighboring sides of the pixel (Fig.fTTI 
b). The derivative dA/dP can be computed with help of the inward directed normal vectors rij at crossing points and 
the vector of the projected passage A (Fig.[TT]c); 



dA njA 

5? "4 Mi 



(14) 



where the index /' runs for all (0, 1 or 2) crossing points belonging to the line segment. For outliers the derivative is 
the unit vector 



3.5. Minimization 

The joint chi-square for all channels in a hit cluster can be written as 



X 



'(£,/') = 24 (A(£,Z,(i'))) 



(16) 



where the index / runs for all the touched and untouched channels of the hit. Note that while the deposits y, are given, 
the As depend on e and the path lengths /, inside the channels (Eq. (|4]i) which in turn depend on the actual hit position 
P. The best hit position, and the most probable differential energy loss will be estimated by minimizing The (three 
dimensional) path lengths /, can be obtained from the projected lengths /(, as 

A 

where / is the total path length and A is the total projected length, both are fixed from the local trajectory direction. 
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Figure 12: Residuals of the reconstructed hit position in directions perpendicular to (left) or parallel with (right) the projected trajectory using the 
discussed methods. 
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Figure 13: Resolution of the reconstructed liit position in directions perpendicular to (left) or parallel with (right) the projected trajectory using the 
discussed methods, as a function of the number of pixels in the cluster. Lines are drawn to guide the eye. 



In case of hit position estimation the value of e is not relevant, hence instead of the complete relationship shown 
in Eq. ^ a simple linearity can be assumed 

A(0 ^el (17) 
where e is rather the average rate of energy loss. The first derivatives of the components with respect to (e, P) are 

dx^ _ dx^^ dx^ _ dx^ ^ dl 

de ~ dA ' dPj ~ dA ^ dPj 

while the second derivatives are 

where the term containing -1^ = is not shown. The second derivative cross-term is 

dsdPj \<9A2 dAjdPj' 

The partial derivative of the path length is 

dl _ I dA 
dP ~ AdP 

and the last term is given in Eqs. (fT4l l and ( fTSl ). During minimization the positivity of s can be assured by choosing 
logs as the free parameter Since both the first and second derivatives can be calculated, Newton's method (l)^ can 
be used for fast and precise minimization. 

In case of strips the observed signals are correlated. Before starting with the position estimation the coupling has 
to be undone first: it will be discussed in Sec. 14. II 

3.6. Results 

For the studies in Sec.[3]and|4]charged particles in the /3y range of 0.56 - 10.0 were generated with flat distribution 
in log(j6y), and in pseudo-rapidity, -2.5 < rj < 2.5, following a pr expi-pr/T) shape for the transverse momentum 
distribution, where T - 0.15 GeV/c was set. Some relevant quantities used in the simulation are given in Tabled 



d^X _ 5V 2 dl 
dPjdPk ~ M^"*" WjdP'k 
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Figure 14: Residuals of the reconstructed hit position (left) and the resolution of the reconstructed hit position using the discussed methods, as a 
function of the number of strips in the cluster (light). Lines are drawn to guide the eye. 



Table 2: Some relevant quantities used in silicon detector simulation. 

Detector thickness 300 /vm 

Pitch, pixels {x direction) 100 jim 

Pitch, pixels {y direction) 200 jim 

Pitch, strips 100 jim 



Signal coupling, strips 





.1 


Channel noise 


1.5 


keV 


Channel threshold 


7.5 


keV 


Channel saturation 


150 


keV 
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Figure 15: Residuals of the reconstructed hit position for pixels (x and y directions) and strips, from left to right. The measured values from 
comparison of reconstructed and simulated positions (Fitter, dash-dotted line) are shown together with the distribution predicted from the errors of 
the fitted parameters (Predicted, dashed line). 



The pixel and strip layers were located in a barrel geometry at 10 cm and 50 cm radial distances from the beam axis, 
respectively, in a B = 4 T magnetic field. 100 000 complete pixel and strip hits were generated, down to the level 
of individual channels. For simplicity, the charge drift direction was assumed to be parallel with the electric field, 
perpendicular to the silicon unit. Hence the change in direction due to the E x B effect was omitted. 

The addition of local particle direction and energy loss information will make the residuals of the hit position 
measurement smaller In case of pixels the deviations from the true value can be decomposed into projections parallel 
and perpendicular to the passage of the particle (Fig.[T2li. In both cases the proposed model (labeled as Fitter) gives 
better results than the weighted or first-last methods introduced in Sec. 13.11 The differences are mostly seen in the 
perpendicular direction. The measured Gaussian resolutions as a function of the number of pixels in the cluster for 
all three discussed methods are shown in Fig. [T3] in some cases reaching 1 fim levels. The corresponding plots of 
residuals and position resolution for strips are shown in Fig. [141 again the fitter giving the best results. 

The resolution of the position estimate can be calculated hit by hit using the observed Fisher information (Sec. 12.4b . 
Residuals of the reconstructed hit position for pixels and strips are shown in Fig. [15] The measured values from 
comparison of reconstructed and simulated positions are plotted together with the distribution predicted from the 
errors of the fitted parameters. There is a very good agreement between observed and predicted values. Note that for 
pixels the cross-correlation term can also be deduced giving valuable input for track refit. 



4. Estimation of hit energy deposit 

The deposited energy Q in the cluster could be plainly estimated with the sum of the individual channel deposits 

i 

but this approach clearly has various biases (Sec. l2.21 l. 

The deposited energy can be lost during clustering, because the deposit can be below the channel threshold (left 
truncation). The average lost deposit {y} can be estimated using the discussed microscopical model. Example values 
as function of path length / for several threshold settings (f = 6, 8 and 10 keV) are shown in Fig. [TSlleft. for particles 
with ySy = 3.16. In case of large I the average lost deposit is close to the threshold, (y) x t, because the probability 
density function p{y\l) is very steep for low y values. If the path length I is small the average lost deposit will be simply 
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Figure 16: The average lost deposit (y> as function of the path length / for several threshold. Left: deposits below threshold with / = 6, 8 and 
10 keV, for particles with py = 3.16. Right: deposits above censoring level with t = 100, 150 and 200 keV, for particles with fiy = 0.32. The curves 
show the functional forms of Eq. (18) and Eq. respectively. 




Figure 17: Residuals of the reconstructed deposited energy of the hit cluster for pixels (left) and strips (right). The result of the simple sum and the 
proposed fitter are compared. 



(y) « eI. These observations can be successfully described and matched with simulated data (Fig. [T6}left) using the 
following formula: 

<3;),<, ^ftanhf^V (18) 



_ t _ 

The measured value can also be saturated, when the deposit is above censoring level (right censoring). Example 
values as function of path length / for several threshold settings {t - 100, 150 and 200 keV) are shown in Fig.fTSlright. 
for particles with fiy - 0.32. In case of small I the average lost deposit can be estimated as 

r3/exp[-vy/o-A(f)] o-A(f) 

iy) ~ —n -t + . 

j exp [-t7/o-a(0] ^ 
15 



If the path length I is large the average lost deposit will be simply (y) » si. These observations can be successfully 
described and matched with simulated data (Fig.fTSlright) using the following formula: 



{y}y- 



>•>' 



max 1 1 H , si 



(19) 



4.1. Coupled signals 

In case of strips the observed signals are correlated due to the capacitive coupling of the neighboring strips and 
cross-talk. This effect can be modelled by a tridiagonal matrix: 



C = 



1 — 2a a 
a I —2a 









a 

I -2a) 



where a > 0, such that the measured deposits y' can be obtained from the original ones y by 

y' = Cy. 

Since the signals are coupled, YjiX^ would sum dependent values making the minimization false. The coupling 
should be first undone by applying C"' to the measured deposits. The inverse of C can be exactly calculated u2m . 
furthermore if or <K 1 it is well approximated by 



C 



(2-1 /a) 



Note that due to the presence of signals below threshold and saturated values the decoupling is not exact. The 
measured deposit in the ith channel is either given, or subject to inequalities: 



(Cy)i = y'l 



(Cy)i < t 



{Cy)i > t. 



In that sense a system of equalities and inequalities has to be solved fory. The problem can be easily handled with 
the tools of linear programming Ulllll, the OOQP package [14], available in ROOT [ISJ under MATH/QUADP, 
could be used for this purpose. The solution is not unique, hence the minimization of a scalar product cy' is also 
required. In this case it is advantageous to set c =(!,..., 1). This way those solution will be selected where 
Z/J/ = Qpittei- is minimal. 

In fact the situation is even more complicated due to the presence of detector and readout noise. Assuming a 
Gaussian noise with standard deviation cr, the contribution of the /th channel is 



{Cy)i 



XfKt 



0, 



if 



2 <0 



otherwise 

if M£ik + 2 > 
otherwise. 



Here, in case of the inequalities, the exact values 

P(y < t\m) 



1 + erf 



■V2 
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were approximated as 



X = -21ogP <^ 2 \ ,T ; ' 

lO, otherwise. 

In order to keep the minimization away from unphysical negative y, values, a penalty (Cj),7cr)2 is added if (Cy)i < 

0. 

The joint chi-square YuX\' ^ positive definite quadratic function of the original deposits y,. The first and second 
derivatives with respect to components are easily calculable and the function can be minimized with Newton's 
method, in few iterative steps. 

If some y, values are negative at the minimum, the one with lowest value is fixed to 0, the minimization is redone. 
In the end we have a minimum joint chi-square with all the y-, values being non-negative. At the minimum the total 
charge and its variance are given as 

respectively. 
4.2. Results 

The details of the simulation are the same as they were described in Sec. 13.61 Residuals of the reconstructed 
deposited energy of the hit cluster are shown in Fig. [17] for pixels and strips. The result of the simple sum and the 
proposed fitter are compared. The sum has a -5 keV shift for both configurations which is a 5% effect for minimum 
ionizing particles along a 300 yum path in silicon. At the same time the fitter has practically no bias and a better 
resolution. 



5. Estimation of the most probable differential energy loss for tracks 

Having the proper hit energy deposits y,, the next step is to estimate the most probable differential energy loss for 
the whole trajectory. The task is very similar (Eq. (fTSIl). but only e has to be optimized, because the path lengths Z, are 
given from track finding. 

Since the terms on the right contain linear or positive definite quadratic functions (Eqs. (fTTTl-(fT3Tl) a fast convergence 
with Newton's method is guaranteed. 
The derivatives are 

dA ifA>y-ycr^(y) 

g;r?<,(A) _ [O, ifA<f-o-A(f) 
dA "l ^'V.w--"" , ifA>f-^A(0 

Xy>X^) _ [-^, ifA<t + crdt) 

dA |o, if A > f -H o-A(t). 

Only the following second derivatives are non-zero: 



5A2 cr^iy)^ 



M2 o-a(0^' 



if A > y - vo-Aiy) 
if A > f - 0-^(0- 
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Figure 18: The expectation value of (t ^ as function of path length / at several fiy values (points). The curves show the functional form defined in 
Eq. )20t . The corresponding fitted powers for e and I are also indicated. 



5.7. Resolution 

The dependence of the standard deviation of the s estimate can be obtained as (Sec. l2.4l i: 

2^ de^ 



where the second derivative is 



1^ 
2 de^ 



\ o-aCv) ) 



if in the Gaussian part 
if in the exponential part. 



According to the fit shown in Fig.[T8]the expectation value of cr^ is a power function of e and I: 



(20) 



With that, also approximating the factor a log(Z,7Zo) in the range / = 50 - 1000 jt/m. 



It is clear that the relative resolution o-(s)/s only slightly depends on e, hence loge is a convenient and uniform 
estimator. Since the exponent of the path length is also small (0.2), it is not the total thickness of the silicon, but the 
number of independent measurements that matters. For a given total path length the relative resolution is proportional 



to n , where n is the number of independent measurements. The reason for that originates in the non-Gaussian 
nature of the energy deposit distribution. Grouped measurements would have worse resolution unlike in the Gaussian 
case where the resolution would stay unchanged. 



5.2. False hit removal 

Since the association of hits to trajectories is not always unambiguous some hits do not belong to the proper track. 
Although their measured deposit is correct, the calculated path length can be false. Assuming that there is at most one 
false hit in a trajectory it can be detected and removed: the exclusion of the outlier hit decreases the joint chi-square 
of the trajectory by a considerable amount. 
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Figure 19: Distribution of tlie estimated most probable differential energy loss values at reference length (e), obtained with the four detector layer 
settings. Several methods are shown: this maximum likelihood fitter, truncated mean, power mean (power -2), harmonic mean (-1) and arithmetic 
mean (1). 



5.3. Distribution of the estimator 

Although the most probable value of e is estimated, it is equally important to deduce, or also estimate, its distribu- 
tion. Although the variance of the estimate is given, the original probability density is a skewed function: the estimator 
is not expected to follow a perfect Gaussian distribution. The distribution of the joint could be determined track by 
track with detailed Monte Carlo simulation, but this possibility would have huge computational demands and timing 
issues. 

Are there other methods that could be used to estimate the distribution of the estimator? One of them is the 

! 1 

Statistical bootstrap [16], that is a random sampling with replacement from the original dataset. The so called jackknife 
method lfl6ll can also be considered where the estimate is systematically recomputed leaving out one observation at a 
time from the sample. Both methods can be used to estimate bias, variance and the shape of the distributions. It can be 
shown that the jackknife method cancels biases with terms proportional to 1 /n and l/n^. A problem common to both 
resampling methods is that they do not work if the sample size is small. Since here we deal sometimes with tracks with 
as few as two or three measurements these approaches are not viable. For example, in case of the jackknife method, 
the variance of the subsamples will have large errors if the subsample consist of only one or two measurements. 
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Figure 20: Left: pion-kaon separation power for four detector layer settings, using several methods: this maximum likelihood fitter, truncated 
mean, power mean (power -2), and harmonic mean (-1). Right: most probable differential energy loss values at reference length (log e), shown for 
the four detector layer settings. The horizontal an'ows show the expected theoretical value. Lines are drawn to guide the eye. 



As a solution, the shape of the log e distribution could be determined by energy deposit regeneration. A measured 
track is used to construct the shape distribution of a given particle species (mass m), if its measured e value is compat- 
ible with the corresponding value at a given pjm. While all kinematical parameters and path lengths of hits in silicon 
are kept, all energy deposit values are randomly regenerated using assumed pjm values, as described in Sec. 12.51 The 
chosen procedure ensures meaningful shape determination even for tracks with very few hits and exploits the success 
of the energy loss model seen at the hit level. 

5.4. Results 

In order to demonstrate track level applications a very simple detector model was used: 16 layers of silicon at 
normal incidence. The orientation of layers and the bending of particle trajectories are taken into account by defining 
four different layer settings: 

• all layers are 300 //m thick (labeled as 300 //m); 

• four groups of layers, each containing four 300, 400, 500 and 600 ^va thick layers (300-600 fim); 

• four groups of layers, each containing four 300, 500, 700 and 900 /im thick layers (300-900 yum); 

• four groups of layers, each containing four 300, 600, 900 and 1200 /im thick layers (300-1200 yum). 

For the study the 100 000 pions and 30 000 kaons were generated at total momentum p - 0.8 GeV/c, and a 
Gaussian noise with 1 keV standard deviation was added to each energy deposit. 

Distribution of the estimated most probable differential energy loss values at reference length (log e), obtained 
with the four detector layer settings, are shown in Fig. [19] (Fitter). Results of several other methods working with the 
differential energy deposit values y,//, are also plotted: truncated mean (average of the lowest half of they,//,- numbers. 



0-50% truncation) lll7Lll8[l : as well as those using all measurements, such as power mean (power -2), harmonic mean 
(power -1) and arithmetic mean (power 1). While the maximum likelihood fitter gives the best results in all the cases, 
power means have much worse resolution. The truncated mean performs surprisingly well. 

The resulted distributions can be fitted with sum of Gaussians. The performance can be further examined by 
plotting the pion-kaon separation power, defined as \m\ - m2\l{{cr\ + cr^ll}^^^, where m,- and cr,- are the means and 
standard deviations for pions and kaons, respectively. Comparisons including estimators listed above are shown in 
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Fig. l20l-left. for the four detector layers settings. Again the maximum likelihood fitter gives the best result over 
the others. The separation power increases if longer hit path lengths or thicker detectors are used (compare with 
observations in Sec. 15. lb . 

Ideally the mean of the log e estimates should not depend on path lengths and detector details. The dependence 
of the mean for the four layers settings is given in Fig.l20lright. The horizontal arrows show the expected theoretical 
value. It is clear that the maximum likelihood fitter provides stable means, while the others show a pronounced 
increase. Although for these latter the dependencies could be compensated, in case of tracks with varying path length 
distribution only the proposed method would perform appropriately. 



6. Detector gain calibration with tracks 

In order to determine the multiplicative gain correction g for a chip, terms such as 

xj{g)^xl,(^) 

should be summed for collected hits and the sum minimized by varying g. The partial derivatives ofx^ are 



With that 



dy 
5y2 



where y = gyj should be substituted. 



2^ JW' ifA<y-v.r.W 

+ 1 2(A-y)o-A(A) 



o-Aiy) 
cTAiyy^ 

' 4vfao-A(A) 
2[o-ACy)+3fc(A-3.)]o-A(A) 



if A > y - vo-Aiy) 



if A < y - vo-Aiy) 
if A > y - VCrA(3'). 



8 dy 



g2 Qyl 



y=gyi 



6.1. Complete gain calibration 

The cross-calibration can be performed in the following steps. 

1. With help of a preliminary gain calibration estimate s for each track. Select pion-like tracks and collect the 
values of Pj, path length and deposit of each hit, and store them for every chip separately. For each chip 
minimize the joint chi-square of all selected hits by varying the gains. The minimization is best performed by a 
golden section search 1.1 9i] first, in order to get near the region of the minimum, followed by a refinement using 
Newton's method. 

2. Using the updated gains select only those tracks which are certainly pions, kaons, protons or decay daughters of 
abundantly produced hadrons (Kg — > n^n^, A — > pn^, A — > pn^), or conversion products (y — > e^e"). Collect 
their hits for every chip separately and minimize again the joint chi-square chip by chip by varying the gains 
with similar methods as above. 

3. Finally, using all tracks, minimize the their joint chi-square simultaneously. This can be accomplished by nested 
minimizations. For a given set of gains the e,- values are determined for each track separately in the course of a 
local minimization. The gains of chips can be highly correlated since tracks often traverse detector units which 
are behind each other or very close in space (double-sided units), but also due to manufacturing details. In 
case of 10"^- 10^ chips the covariance matrix needed for the minimization step is sparse, but huge. It is usually 
numerically difficult to invert. In order to have a treatable problem the detector units with highest correlation 
are identified, this way forming a covariance matrix with a block diagonal structure. This latter is already easily 
inverted and the Newtonian step can be computed. 
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Figure 21: Reconstructed relative gains as function of the real values, shown for the 16 layers used in the detector simulation. 



6.2. Results 

The details of the simulation are the same as they were described in Sec. 15.41 but now the initial gains of the layers 
were set randomly in the range 0.8 - 1 .2. Only the first step in the list of Sec. l6.1l was performed. The extracted relative 
gains as function of the real values are shown in Fig. [21] for all the 16 layers used. While the gains are steadily smaller 
by about 1 % (bias), they have an excellent relative precision. 



7. Conclusions 

In this work a new analytical energy loss model for charged particles in siUcon was introduced. It has few param- 
eters and it is based on the recognition of a special connection between the distribution of the deposited energy and 
the most probable energy loss. Its use was demonstrated through several examples. With help of measured charge 
deposits in pixels or strips of hit clusters their position and energy can be estimated with better accuracy and much 
less bias. Deposits below threshold and saturated values are treated properly, resulting in a wider dynamic range. The 
model was successfully applied to track differential energy loss estimation and to detector gain calibration tasks, again 
showing a performance superior to standard methods. 
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